/* Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

/*
  <file> ieee754.c </file>
  <brief>
    An implementation modf() for 'long double'.
    This file is part of ULS, Unified Lexical Scheme.
  </brief>
  <author>
    Stanley J. Hong <babel96@gmail.com>, May 2011.
  </author>
*/

#include "common.h"
#include "ieee754.h"

#define DOUBLE_SIZE_BYTES     8
#define DOUBLE_EXPO_SIZE_BITS  11
#define DOUBLE_EXPO_BIAS      1023

#define DOUBLE_SIZE_BITS   (DOUBLE_SIZE_BYTES*8)
#define DOUBLE_MENTISA_STARTBIT  (1+DOUBLE_EXPO_SIZE_BITS)

#if LDOUBLE_MEMSIZE == 8
#define LDOUBLE_SIZE_BYTES     DOUBLE_SIZE_BYTES
#define LDOUBLE_EXPOSIZE_BITS  DOUBLE_EXPO_SIZE_BITS
#define LDOUBLE_EXPO_BIAS      DOUBLE_EXPO_BIAS
#elif LDOUBLE_MEMSIZE == 12
#define LDOUBLE_SIZE_BYTES     10
#define LDOUBLE_EXPOSIZE_BITS  15
#define LDOUBLE_EXPO_BIAS      16383
#elif LDOUBLE_MEMSIZE == 16
#define LDOUBLE_SIZE_BYTES     16
#define LDOUBLE_EXPOSIZE_BITS  15
#define LDOUBLE_EXPO_BIAS      16383
#else
#error "unsupported size of long double!"
#endif

#define LDOUBLE_SIZE_BITS   (LDOUBLE_SIZE_BYTES*8)
#define LDOUBLE_MENTISA_STARTBIT  (1+LDOUBLE_EXPOSIZE_BITS)

void
ieee754_bankup_bits(uls_byte_t* srcptr, int start_bit, int end_bit, int n_shift)
{
	uls_uint32  i1, i2, i1_dash, i2_dash;
	uls_uint32  j1, j2,  j1_dash, j2_dash, j;
	uls_uint32  l1, l2, n_bits, n_bytes, m, k;

	uls_byte_t ch_src, ch_carry, mask, mask_carry;
	uls_byte_t* dstptr, *srcptr0;

	//
	//  This shifts left by 'n_shift' the bits of range [start_bit,end_bit] from 'srcptr'.
	//  assert: (end_bit + 1) % 8 == 0
	//
	//     A0         A0+1        A0+2        A0+3        A0+4
	//     ------------------------------------------------------------------
	//     |        ###|$$=========|           |      ###$$|**********|
	//     --------|---|--|--------------------|-----|-----|-----------------
	//             P2  P2"P2'                        P1    P1'
	//             i2     i2'                        i1    i1'
	//             j2     j2'                        j1    j1'
	//
	//   P1(i1,j1) == P(start_bit/8,start_bit%8) :: P1'(i1',j1')
	//   P2(i2,j2) :: P2'(i2',j2')
	//
	//   Let length(P1-P1') == l1 (bits) and
	//     length(P2"-P2) == l2 (bits)
	//

	srcptr0 = srcptr;

	uls_assert((end_bit + 1) % 8 == 0);
	if (n_shift == 0 || start_bit > end_bit) return;
	if (start_bit < n_shift) n_shift = start_bit;

	i1 = start_bit / 8;
	j1 = start_bit % 8;
	m = start_bit - n_shift;
	i2 = m / 8;
	j2 = m % 8;

	n_bits = end_bit - start_bit + 1;
	n_bytes = n_bits / 8;

	if (j1 > 0) {
		l1 = 8 - j1;
		if (j2 > 0) {
			l2  = 8 - j2;
			if (l2 > l1) {
					i2_dash = i2;
					j2_dash = j2 + l1;
			} else { // l2 <= l1
				i2_dash = i2 + 1;
				j2_dash = l1 - l2;
			}
		} else  {
			i2_dash = i2;
			j2_dash = l1;
		}

		i1_dash = i1 + 1;
		j1_dash = 0;

		dstptr = srcptr + i1;
		for (j=j1; j<8; j++) {
			if (uls_test_bit(j, dstptr)) {
				uls_set_bit(m+j-j1, srcptr);
			} else {
				uls_clear_bit(m+j-j1, srcptr);
			}
		}
	} else {
		i2_dash = i2;
		j2_dash = j2;
		i1_dash = i1;
		j1_dash = 0;
	}

	if (n_bytes > 0) {
		dstptr = srcptr + i2_dash;

		srcptr += i1_dash;

		if (j2_dash > 0) {
			// copy the bytes left
			//                 (i2',i2'+1  <--  (i1')
			//               (i2'+1,i2'+2) <--  (i1'+1)
			//               (i2'+2,i2'+3) <--  (i1'+2)
			//                  ........   <--  .......
			// (i2'+n_bytes-1,i2'+n_bytes) <-- (i1'+n_bytes-1)

			m = 8 - j2_dash;
			mask =  ((1 << m) - 1) << j2_dash;
			mask_carry = ((1 << j2_dash) - 1) << m;
			ch_carry = dstptr[0] & mask_carry;

			for (k=0; k<n_bytes; k++) {
				ch_src = srcptr[k];
				dstptr[k] = ch_carry | ((ch_src & mask) >> j2_dash);
				ch_carry = ch_src << m;
			}

			dstptr[k] = (dstptr[k] & ~mask_carry) | ch_carry;

		} else {
			memmove(dstptr, srcptr, n_bytes);
		}
	}

	// and clear the bits trailing the shifted bits.
	uls_clear_bits(srcptr0, end_bit - n_shift + 1, end_bit);
}

static int
get_ieee754_biased_expo(uls_byte_t* buff, int n_bits_expo)
{
	uls_byte_t  mask, ch, ch_carry;
	int expo, i, j1, j2, n_bits, m;

	ch_carry = buff[0] & 0x7F;

	for (expo=0, i=0; i<n_bits_expo/8; ) {
		++i;
		ch = ch_carry << 1;
		if (buff[i] & 0x80) {
			ch |= 0x01;
		}

		expo = (expo << 8) | ch;
		ch_carry = buff[i] & 0x7F;
	}

	if ((n_bits = n_bits_expo % 8) > 0) {
		j1 = 1;
		j2 = n_bits;
		m = 7 - j2;

		mask = ((1 << n_bits) - 1) << m;
		ch = (buff[i] & mask) >> m;
		expo = (expo << n_bits) | ch;

	} else {
		expo = (expo << 7) | ch_carry;
	}

	return expo;
}

static void
put_ieee754_biased_expo(int m_expo, uls_byte_t* buff, int n_bits_expo)
{
	uls_byte_t  mask, ch;
	int i, n_bits, m;

	i = n_bits_expo / 8;

	if ((n_bits = (1+n_bits_expo) % 8) > 0) {
		m = 8 - n_bits;
		mask = (1 << n_bits) - 1;

		ch = (m_expo & mask) << m;
		mask <<= m;

		buff[i] = (buff[i] & ~mask) | ch;
		--i;
		m_expo >>= n_bits;
	}

	for ( ; i>0; i--, m_expo>>=8) {
		buff[i] = m_expo & 0xFF;
	}

	mask = 0x7F;
	buff[0] = (buff[0] & ~mask) | (m_expo & mask);
}

int
ieee754_longdouble_isspecial(long double x, char* nambuf)
{
	// ret-val == 1 : x is a special-value with nambuf filled.
	// ret-val == 0 x is a finite value
	uls_byte_t  buff[LDOUBLE_SIZE_BYTES];
	int expo_biasd, stat = 0;
	int i, i0, minus;

	memcpy(buff, &x, LDOUBLE_SIZE_BYTES);
	arch2be_array(buff, LDOUBLE_SIZE_BYTES);

	minus = (buff[0] & 0x80) ? 1 : 0;

	expo_biasd = get_ieee754_biased_expo(buff, LDOUBLE_EXPOSIZE_BITS);
	i0 = LDOUBLE_MENTISA_STARTBIT;

	if (expo_biasd == (1<<LDOUBLE_EXPOSIZE_BITS)-1) {
		if (uls_find_first_1bit(buff,
				i0, LDOUBLE_SIZE_BITS-1, (uls_uint32 *) &i) == 0) {
			if (minus) strcpy(nambuf, "-INF");
			else strcpy(nambuf, "+INF");

		} else if (i == i0) {
			strcpy(nambuf, "QNaN");
		} else { // i > i0
			strcpy(nambuf, "SNaN");
		}
		stat = 1;

	} else if (expo_biasd == 0) {
		if (uls_find_first_1bit(buff,
			i0, LDOUBLE_SIZE_BITS-1, (uls_uint32 *) &i) == 0) {
			if (minus) strcpy(nambuf, "-0");
			else strcpy(nambuf, "+0");
		} else {
			strcpy(nambuf, "PDR");
		}
	} else {
		strcpy(nambuf, "PNR");
	}

	return stat;
}

long double
ieee754_modlf(long double x, long double* p_frac)
{
	uls_byte_t  buff[LDOUBLE_SIZE_BYTES];
	uls_byte_t  buff2[LDOUBLE_SIZE_BYTES];
	int m_expo, i0, i;

	memcpy(buff, &x, LDOUBLE_SIZE_BYTES);
	arch2be_array(buff, LDOUBLE_SIZE_BYTES);

	m_expo = get_ieee754_biased_expo(buff, LDOUBLE_EXPOSIZE_BITS) - LDOUBLE_EXPO_BIAS;
	if (m_expo < 0) {
		if (p_frac!=NULL) *p_frac = x;
		return 0.;
	}

	// m_expo >= 0
	i0 = LDOUBLE_MENTISA_STARTBIT + m_expo;

#ifdef ULS_IMPLICIT_1BIT_MENTISA
	if (p_frac != NULL) {
		memcpy(buff2, buff, LDOUBLE_SIZE_BYTES);

		if (uls_find_first_1bit(buff2, i0, LDOUBLE_SIZE_BITS-1, (uls_uint32 *) &i)==0) {
			*p_frac = 0.;
		} else {
			if (++i >= LDOUBLE_SIZE_BITS) {
				uls_clear_bits(buff2, LDOUBLE_MENTISA_STARTBIT, LDOUBLE_SIZE_BITS-1);
			} else {
				ieee754_bankup_bits(buff2, i, LDOUBLE_SIZE_BITS-1, i-LDOUBLE_MENTISA_STARTBIT);
			}
			put_ieee754_biased_expo(i0-i+LDOUBLE_EXPO_BIAS, buff2, LDOUBLE_EXPOSIZE_BITS);
			arch2be_array(buff2, LDOUBLE_SIZE_BYTES);
			memcpy(p_frac, buff2, LDOUBLE_SIZE_BYTES);
		}
	}

	uls_clear_bits(buff, i0, LDOUBLE_SIZE_BITS-1);
#else
	if (p_frac != NULL) {
		memcpy(buff2, buff, LDOUBLE_SIZE_BYTES);

		if (uls_find_first_1bit(buff2, i0+1, LDOUBLE_SIZE_BITS-1, (uls_uint32 *) &i)==0) {
			*p_frac = 0.;
		} else {
			ieee754_bankup_bits(buff2, i, LDOUBLE_SIZE_BITS-1, i-LDOUBLE_MENTISA_STARTBIT);
			put_ieee754_biased_expo(i0-i+LDOUBLE_EXPO_BIAS, buff2, LDOUBLE_EXPOSIZE_BITS);
			arch2be_array(buff2, LDOUBLE_SIZE_BYTES);
			memcpy(p_frac, buff2, LDOUBLE_SIZE_BYTES);
		}
	}

	uls_clear_bits(buff, i0+1, LDOUBLE_SIZE_BITS-1);
#endif

	arch2be_array(buff, LDOUBLE_SIZE_BYTES);
	memcpy(&x, buff, LDOUBLE_SIZE_BYTES);

	return x;
}
